Chromosome compaction - simple three-state model

library(ChromCom)
library(ggplot2)
library(gridExtra)
library(reshape2)

11/07/2017

Writing introduction. First script to do simulation.

12/07/2017

Let’s try if it works. An example for \(t_1 = -60\) min, \(\Delta t = 0\), \(k_1 = 0.04\ {\rm min}^{-1}\) and \(k_2 = 0.04\ {\rm min}^{-1}\).

pars <- list(
  t1 = -60,
  dt = 0,
  r1 = 0.04,
  r2 = 0.04
)
chr <- ChromCom3(pars)
chr <- generateCells(chr, nsim=1000)
plotTimelines(chr)

Just to check if output is as expected. With \(k_2=0\) we should have pure exponential decay. The yellow curve is \(e^{-k_1 t}\).

pars <- list(
  t1 = -60,
  dt = 0,
  r1 = 0.04,
  r2 = 0
)
tchr <- ChromCom3(pars)
tchr <- generateCells(tchr, nsim=1000)
x <- seq(from=pars$t1, to=100, by=1)
y <- exp(-pars$r1 * (x - pars$t1))
g <- plotTimelines(tchr)
g + geom_line(data=data.frame(x=x, y=y), aes(x,y), colour="yellow")

Close, but not perfect. Here is an example with much smaller time step (0.1) and larger number of cells (10,000).

pars <- list(
  t1 = -60,
  dt = 0,
  r1 = 0.04,
  r2 = 0
)
tchr <- ChromCom3(pars, time = seq(from=-140, to=90, by=0.1))
tchr <- generateCells(tchr, nsim=10000)
x <- seq(from=pars$t1, to=100, by=1)
y <- exp(-pars$r1 * (x - pars$t1))
g <- plotTimelines(tchr)
g + geom_line(data=data.frame(x=x, y=y), aes(x,y), colour="yellow")

Parameter grid

parameterGrid <- function(pars, par1, range1, par2, range2) {
  melts <- NULL
  for(p1 in range1) {
    pars[[par1]] <- p1
    for(p2 in range2) {
      pars[[par2]] <- p2
      chr <- ChromCom3(pars)
      chr <- generateCells(chr)
      label1 <- sprintf("%s = %.3g", par1, p1)
      label2 <- sprintf("%s = %.3g", par2, p2)
      m <- meltTimelines(chr, label1=label1, label2=label2)
      melts <- rbind(melts, m)
    }
  }
  timelinePanel(melts)
}

pars <- list(
  t1 = -30,
  dt = 0,
  r1 = 0.05,
  r2 = 0.01
)
range1 <- c(0.01, 0.05, 0.08, 0.12)
range2 <- c(0.01, 0.05, 0.08, 0.12)
parameterGrid(pars, "r1", range1, "r2", range2)

I think we need to introduce a delay between pink and red.

pars <- list(
  t1 = -30,
  dt = 10,
  r1 = 0.05,
  r2 = 0.01
)
range1 <- c(0.01, 0.03, 0.05, 0.10)
range2 <- c(0.01, 0.05, 0.08, 0.12)
parameterGrid(pars, "r1", range1, "r2", range2)

What about range of delays?

pars <- list(
  t1 = -30,
  dt = 0,
  r1 = 0.03,
  r2 = 0.01
)
range1 <- c(0, 10, 20, 30)
range2 <- c(0.01, 0.05, 0.08, 0.12)
parameterGrid(pars, "dt", range1, "r2", range2)

Experimental data.

echr <- experimentalData(dataFile$scramble)
plotTimelines(echr, smooth=TRUE, k=15)

13/07/2017

Prepared shiny app for deployment.

24/07/2017

Tomo suggested that the model parameters can be derived directly from data. We can write the following equations

\(\dot{B} = -k_1 B\)

\(\dot{P} = k_1 B - k_2 R\)

\(\dot{R} = k_2 P\)

\(B + P + R = 1\)

However, I’m not sure if this is exacly the model I use in the simulation. In the simulation, we assumed that P->R transition can happen only after the BB->P transition.

Because proportions and their derivatives are observed, we can find the rates:

\(k_1 = - \frac{\dot{B}}{B}\)

\(k_2 = \frac{\dot{R}}{P}\)

This will require a lot of smoothing, otherwise the derivatives will be all over the place.

First, I see if I can recover rates from the generated data.

timeDeriv <- function(chr, k=20) {
  ts <- list()
  for(col in chr$colours) {
    x <- chr$cnt[[col]] / chr$cnt$total
    x <- caTools::runmean(x, k)
    x <- ts(x, start=chr$timepars$start, deltat=chr$timepars$step)
    xdot <- diff(x)
    ts[[col]] <- x
    ts[[paste0(col, ".diff")]] <- xdot
  }
  ts
}

plotRates <- function(chr, k=20) {
  ts <- timeDeriv(chr, k=k)
  k1 <- -ts$BB.diff / ts$BB
  k2 <- ts$R.diff / ts$P
  k1[which(is.nan(k1))] <- NA
  k2[which(is.nan(k2))] <- NA
  df <- data.frame(t=time(k1), k1=k1, k2=k2)
  m <- melt(df, measure.vars = c("k1", "k2"))
  m$value <- as.numeric(m$value)
  m$value[which(m$value > 0.5)] <- 0.5
  m$value[which(m$value < -0.5)] <- -0.5
  ggplot(m, aes(x=t, y=value)) + geom_line(aes(colour=variable)) 
}

Start with a model wiht 10,000 cells and 0.1 min time step. In this model \(k_1 = 0.04\) and \(k_2 = 0\). Top panels show original data, lower panels - smoothed data.

grid.arrange(plotTimelines(tchr), plotRates(tchr, k=1), ncol=2)

grid.arrange(plotTimelines(tchr, smooth=TRUE, k=200), plotRates(tchr, k=200), ncol=2)

What happens if I use only 100 cells and 1 min time step in the model? This time \(k_2 = 0.04\).

chr <- generateCells(chr, nsim=100)
grid.arrange(plotTimelines(chr, k=1), plotRates(chr, k=1), ncol=2)

grid.arrange(plotTimelines(chr, smooth=TRUE, k=20), plotRates(chr, k=20), ncol=2)

And now, for real data. Here are scramble data smoothed with running mean with window size of 10 time points.

grid.arrange(plotTimelines(echr, smooth=TRUE, k=10), plotRates(echr, k=10), ncol=2)

Now, smoothing with 80 points.

grid.arrange(plotTimelines(echr, smooth=TRUE, k=80), plotRates(echr, k=80), ncol=2)